#########################################
# Replication file for Neuquen Province #
#########################################


rm(list=ls(all=TRUE))
#setwd("")
library(eiPack)
library(coda)
library(foreign)

data=read.csv(file.path(getwd(),"data","NQN_data.csv"))


# Define row and column marginals of the transition matrix. 
GMPN=data$gral_153
GCCN=data$gral_532
GFPV=data$gral_501
GO=data$gral_505+data$gral_40+data$gral_533+data$gral_166+data$gral_23
GN=data$votantes-GMPN-GCCN-GFPV-GO


P_AMPN=data$paso_3073
P_BMPN=data$paso_3072
PCCNA=data$paso_3077
PCCNB=data$paso_3076
PFPV=data$paso_501
PO=data$paso_505+data$paso_40+data$paso_533+data$paso_166+data$paso_23+data$paso_183
PNV=data$votantes-P_AMPN-P_BMPN-PCCNA-PCCNB-PFPV-PO

data1=data.frame(P_AMPN,P_BMPN,PCCNA,PCCNB,PFPV,PO,PNV,GMPN,GCCN,GFPV,GO,GN)


# Tunning the EI algorithm
tune.nocov=tuneMD(cbind(GMPN,GCCN,GFPV,GO,GN)~cbind(P_AMPN,P_BMPN,PCCNA,PCCNB,PFPV,PO,PNV),
                  data=data1,ntunes=1,totaldraws=200000)


# MCMC
out.nocov=ei.MD.bayes(cbind(GMPN,GCCN,GFPV,GO,GN)~cbind(P_AMPN,P_BMPN,PCCNA,PCCNB,PFPV,PO,PNV),
                      covariate=NULL,data=data1,tune.list=tune.nocov,verbose=0,ret.beta="r",thin=80,burnin=100000,ret.mcmc=TRUE) 




# Obtain transition matrix for each mesa: 
source(file.path(getwd(),"Codes","getbetas.R"))
beta11=getbetas(out.nocov,"P_AMPN","GMPN")
beta12=getbetas(out.nocov,"P_AMPN","GCCN")
beta13=getbetas(out.nocov,"P_AMPN","GFPV")
beta14=getbetas(out.nocov,"P_AMPN","GO")
beta15=getbetas(out.nocov,"P_AMPN","GN")
beta21=getbetas(out.nocov,"P_BMPN","GMPN")
beta22=getbetas(out.nocov,"P_BMPN","GCCN")
beta23=getbetas(out.nocov,"P_BMPN","GFPV")
beta24=getbetas(out.nocov,"P_BMPN","GO")
beta25=getbetas(out.nocov,"P_BMPN","GN")
beta31=getbetas(out.nocov,"PCCNA","GMPN")
beta32=getbetas(out.nocov,"PCCNA","GCCN")
beta33=getbetas(out.nocov,"PCCNA","GFPV")
beta34=getbetas(out.nocov,"PCCNA","GO")
beta35=getbetas(out.nocov,"PCCNA","GN")
beta41=getbetas(out.nocov,"PCCNB","GMPN")
beta42=getbetas(out.nocov,"PCCNB","GCCN")
beta43=getbetas(out.nocov,"PCCNB","GFPV")
beta44=getbetas(out.nocov,"PCCNB","GO")
beta45=getbetas(out.nocov,"PCCNB","GN")
beta51=getbetas(out.nocov,"PFPV","GMPN")
beta52=getbetas(out.nocov,"PFPV","GCCN")
beta53=getbetas(out.nocov,"PFPV","GFPV")
beta54=getbetas(out.nocov,"PFPV","GO")
beta55=getbetas(out.nocov,"PFPV","GN")
beta61=getbetas(out.nocov,"PO","GMPN")
beta62=getbetas(out.nocov,"PO","GCCN")
beta63=getbetas(out.nocov,"PO","GFPV")
beta64=getbetas(out.nocov,"PO","GO")
beta65=getbetas(out.nocov,"PO","GN")
beta71=getbetas(out.nocov,"PNV","GMPN")
beta72=getbetas(out.nocov,"PNV","GCCN")
beta73=getbetas(out.nocov,"PNV","GFPV")
beta74=getbetas(out.nocov,"PNV","GO")
beta75=getbetas(out.nocov,"PNV","GN")
betaout=data.frame(beta11,beta12,beta13,beta14,beta15,
                   beta21,beta22,beta23,beta24,beta25,
                   beta31,beta32,beta33,beta34,beta35,
                   beta41,beta42,beta43,beta44,beta45,
                   beta51,beta52,beta53,beta54,beta55,
                   beta61,beta62,beta63,beta64,beta65,
                   beta71,beta72,beta73,beta74,beta75,
                   data$votantes,data$prov,data$comuna,data$circ,data$mesa)
write.table(betaout,file=file.path(getwd(),"Output","NQN_Mesa.txt"),col.names=NA)


# Obtain aggregate transition matrix. 
source(file.path(getwd(),"Codes","getbetasagr.R"))
beta11=getbetasagr(out.nocov,"P_AMPN","GMPN")
beta12=getbetasagr(out.nocov,"P_AMPN","GCCN")
beta13=getbetasagr(out.nocov,"P_AMPN","GFPV")
beta14=getbetasagr(out.nocov,"P_AMPN","GO")
beta15=getbetasagr(out.nocov,"P_AMPN","GN")
beta21=getbetasagr(out.nocov,"P_BMPN","GMPN")
beta22=getbetasagr(out.nocov,"P_BMPN","GCCN")
beta23=getbetasagr(out.nocov,"P_BMPN","GFPV")
beta24=getbetasagr(out.nocov,"P_BMPN","GO")
beta25=getbetasagr(out.nocov,"P_BMPN","GN")
beta31=getbetasagr(out.nocov,"PCCNA","GMPN")
beta32=getbetasagr(out.nocov,"PCCNA","GCCN")
beta33=getbetasagr(out.nocov,"PCCNA","GFPV")
beta34=getbetasagr(out.nocov,"PCCNA","GO")
beta35=getbetasagr(out.nocov,"PCCNA","GN")
beta41=getbetasagr(out.nocov,"PCCNB","GMPN")
beta42=getbetasagr(out.nocov,"PCCNB","GCCN")
beta43=getbetasagr(out.nocov,"PCCNB","GFPV")
beta44=getbetasagr(out.nocov,"PCCNB","GO")
beta45=getbetasagr(out.nocov,"PCCNB","GN")
beta51=getbetasagr(out.nocov,"PFPV","GMPN")
beta52=getbetasagr(out.nocov,"PFPV","GCCN")
beta53=getbetasagr(out.nocov,"PFPV","GFPV")
beta54=getbetasagr(out.nocov,"PFPV","GO")
beta55=getbetasagr(out.nocov,"PFPV","GN")
beta61=getbetasagr(out.nocov,"PO","GMPN")
beta62=getbetasagr(out.nocov,"PO","GCCN")
beta63=getbetasagr(out.nocov,"PO","GFPV")
beta64=getbetasagr(out.nocov,"PO","GO")
beta65=getbetasagr(out.nocov,"PO","GN")
beta71=getbetasagr(out.nocov,"PNV","GMPN")
beta72=getbetasagr(out.nocov,"PNV","GCCN")
beta73=getbetasagr(out.nocov,"PNV","GFPV")
beta74=getbetasagr(out.nocov,"PNV","GO")
beta75=getbetasagr(out.nocov,"PNV","GN")
betaoutagr=data.frame(beta11,beta12,beta13,beta14,beta15,
                   beta21,beta22,beta23,beta24,beta25,
                   beta31,beta32,beta33,beta34,beta35,
                   beta41,beta42,beta43,beta44,beta45,
                   beta51,beta52,beta53,beta54,beta55,
                   beta61,beta62,beta63,beta64,beta65,
                   beta71,beta72,beta73,beta74,beta75,
                   data$votantes,data$prov,data$comuna,data$circ,data$mesa)

betaoutagr=as.vector(betaoutagr[1,1:140])
betaoutagr=matrix(betaoutagr,4)
betaoutagr_coefs=t(matrix(betaoutagr[1,],5))
betaoutagr_sd=t(matrix(betaoutagr[2,],5))
write.table(betaoutagr_coefs,file=file.path(getwd(),"Output","NQNAggregate_Coefs.txt"),col.names=NA)
write.table(betaoutagr_sd,file=file.path(getwd(),"Output","NQNAggregate_SD.txt"),col.names=NA)



# save(out.nocov, file = "NQNMCMC.Rdata") # MCMC runs

# Geweke Convergence Statistics. 
geweke=geweke.diag(out.nocov$draws$Beta)
gewekestats=geweke[[1]]
gewekedf=data.frame(gewekestats)
#write.table(gewekedf,file="NQNgewekelist.txt",col.names=NA)

jpeg(file.path(getwd(),"Output","NQNHistogramGeweke.jpeg"))
hist(gewekestats)
dev.off()

print("Max Geweke")
print(max(gewekestats))
print("Min Geweke")
print(min(gewekestats))
print("Percentage Above 2.35")
print(length(gewekestats[abs(gewekestats)>2.35])/length(gewekestats))
print("Percentage Above 1.96")
print(length(gewekestats[abs(gewekestats)>1.96])/length(gewekestats))
print("Percentage Above 1.5")
print(length(gewekestats[abs(gewekestats)>1.5])/length(gewekestats))
print("Percentage Above 1.5")
print(length(gewekestats[abs(gewekestats)>1])/length(gewekestats))
print("Percentage Above 0.5")
print(length(gewekestats[abs(gewekestats)>0.5])/length(gewekestats))

